clear all

% Coded by Jeronimo Callejas
% First version: May 25, 2021.
% This  version: June 2, 2021.`

% Load data and create variables.

[income,txt2,raw2]=xlsread('inc.csv');

T = readtable('Colombia_final.xlsx','ReadVariableNames',true);
T=sortrows(T,[183,38]);
Z=T(:,48:71); % Set of instruments
Y_M=table2array(T(:,184));


% clear dataset
power=table2array(T(:,196));
consumption=table2array(T(:,195));
high_end=table2array(T(:,197));
num_mo_ft=table2array(T(:,75)); % number of models each firm whithin fuel type
p_nt=table2array(T(:,9));

% remove nan before going to estimation

Temp=[table2array(Z), power, consumption, high_end,num_mo_ft, p_nt];
[rowN, colN] = find(isnan(Temp));
remove=unique(rowN);

T(remove,:)=[];

%%
Z=T(:,48:71); % Set of instruments
Y_M=table2array(T(:,184));
power=table2array(T(:,196));
consumption=table2array(T(:,195));
high_end=table2array(T(:,197));
num_mo_ft=table2array(T(:,75)); % number of models each firm whithin fuel type
p_nt=table2array(T(:,9));
emissions=table2array(T(:,20));
energy_con=table2array(T(:,21));
size_car=table2array(T(:,3));

% Commodities
cobalt=table2array(T(:,201));
lithium=table2array(T(:,202));
steel=table2array(T(:,203));

% Trasporte

trans=table2array(T(:,204));

% Getting dummies

d_make=table2array(T(:,76:135));
d_sta_sg=table2array(T(:,136:143));
d_fue_ty=table2array(T(:,144:154));
d_len_cl=table2array(T(:,155:160));
d_bod_ty=table2array(T(:,161:169));
d_dri_wh=table2array(T(:,170:171));
d_year=table2array(T(:,183));
d_year_month=dummyvar(categorical(Y_M));
d_year_month(:,1)=[];
d_size=dummyvar(categorical(size_car));

% Getting estimation variables

transmission=table2array(T(:,18));
q=table2array(T(:,41));

year=table2array(T(:,183));
tm=table2array(T(:,188));
fuel_ef=table2array(T(:,194));
turbo=table2array(T(:,15));
hp=table2array(T(:,16));
kw=table2array(T(:,200));
ccm=table2array(T(:,32));
weight=table2array(T(:,4));
bat_size=table2array(T(:,198));
drv_rng=table2array(T(:,199));
num_gears=table2array(T(:,33));
num_mo_fs=table2array(T(:,72)); % number of models each firm whithin std seg
num_mo=table2array(T(:,73)); % number of vesrsion same models each firm
num_mo_f=table2array(T(:,74)); % number of models each firm
char_est=table2array(T(:,186)); % number of charging stations
month=table2array(T(:,38));
power2=log(kw./weight);

ss_car=zeros(size(p_nt));
ml_car=zeros(size(p_nt));
mm_car=zeros(size(p_nt));
mb_car=zeros(size(p_nt));
ss_car(ccm<1200)=1;
ml_car(ccm>1201 & ccm<1600)=1;
mm_car(ccm>1601 & ccm<2000)=1;
mb_car(ccm>2001 & ccm<2400)=1;
engine_size2=[ss_car,ml_car,mm_car,mb_car];

small_car=zeros(size(p_nt));
medium_car=zeros(size(p_nt));
small_car(ccm<1601)=1;
medium_car(ccm>1600 & ccm<2801)=1;
engine_size=[small_car,medium_car];

% Create index
N=length(p_nt);
id=categorical(table2array(T(:,1)));
id_ver=categorical(table2array(T(:,19)));
id_make=categorical(table2array(T(:,177)));
id_market=categorical(Y_M);
id_fueltype=categorical(table2array(T(:,179)));
id_fuelcat=categorical(table2array(T(:,37)));
id_staseg=categorical(table2array(T(:,178)));

d_fuelcat=dummyvar(id_fuelcat);
d_fuelcat(:,1)=[];

d_month=dummyvar(month);
d_month(:,1)=[];

d_staseg=dummyvar(id_staseg);
d_staseg(:,1)=[];

d_num_gears=dummyvar(num_gears);
d_num_gears(:,1)=[];

ind_market=grp2idx(id_market);
ind_make=grp2idx(id_make);
ind_ver=grp2idx(id_ver);
QI=dummyvar(id_market);

temp=repmat(ind_make,1,N);
A1=eq(temp,temp');% ownership matrix
clear temp

tariffs=ones(N,3);
tariffs=tariffs.*[0.35,0.05,0];

electric=sum(table2array(T(:,144:145)),2);
hybrid=sum(table2array(T(:,[146:150,152])),2);
non_elec=ones(size(N,1))-electric-hybrid;

elec_char_est=[electric, electric.*char_est];

tariffs=sum(tariffs.*[non_elec, hybrid, electric],2);
tariffs(year<2017)=0.35;

sales_tax=tariffs;
sales_tax(sales_tax~=0.35)=0.05;
sales_tax(sales_tax==0.35)=0.19;

p=(p_nt.*(1+sales_tax))./1000;

edges=[min(power),quantile(power,0.25), quantile(power,0.50), quantile(power,0.75),max(power)];
d_power=dummyvar(discretize(power,edges));
edges=[min(consumption),quantile(consumption,0.25), quantile(consumption,0.50), quantile(consumption,0.75),max(consumption)];
d_consumption=dummyvar(discretize(consumption,edges));

edges_p=[min(p),quantile(p,0.20), quantile(p,0.40), quantile(p,0.60), quantile(p,0.80),max(p)];
d_price=dummyvar(discretize(p,edges_p));

cons_cat=sum(d_consumption.*[1,2,3,4],2);
cons_cat(electric==1)=0;

pr_cat=sum(d_price.*[1,2,3,4,5],2);
%p=p_nt;
%% Create additional instruments
% Competitiors characteristics standard deviation
new_inst_char=[ccm, hp, weight, consumption];
[xx, yy] = ndgrid(ind_market,1:size(new_inst_char,2));
new_inst_sum=accumarray([xx(:) yy(:)],new_inst_char(:));
num_prod=sum(QI)';
num_prod=num_prod(ind_market);
new_inst_mean=(new_inst_sum(ind_market,:)-new_inst_char)./(num_prod-1);
temp_sd=(new_inst_char-new_inst_mean).^2;
new_inst_sd_sum=accumarray([xx(:) yy(:)],temp_sd(:));
new_inst_sd=((new_inst_sd_sum(ind_market,:)-temp_sd)./(num_prod-1)).^0.5;

% Competitiors characteristics standard deviation within standart segment group.
[~, ~, jj] = unique([id_market,id_staseg], 'rows');
[xx, yy] = ndgrid( jj,1:size(new_inst_char,2));
new_inst_sum_seg=accumarray([xx(:) yy(:)],new_inst_char(:));
new_inst_sum_deno=accumarray(jj,ones(size(p)));

new_inst_sum_seg=new_inst_sum_seg(jj,:)-new_inst_char;
new_inst_sum_deno=repmat(new_inst_sum_deno(jj),1,4);

new_inst_sum_mean=new_inst_sum_seg./(new_inst_sum_deno);
temp_sd_seg=(new_inst_char-new_inst_sum_mean).^2;
new_inst_seg_sd_sum=accumarray([xx(:) yy(:)],temp_sd_seg(:));
new_inst_seg_sd=((new_inst_seg_sd_sum(jj,:)-temp_sd_seg)./(new_inst_sum_deno)).^0.5;

% Interaction automatic standard segment
int_tra_seg=transmission.*d_staseg;
int_he_fueltype=high_end.*d_fuelcat;

% Supply side instrumens
bat_x_lit=bat_size.*lithium;
bat_x_cob=bat_size.*cobalt;
wei_x_ste=weight.*steel;
tra_x_wei=trans.*weight;
save('data')
%% Logit - Nested logit
load data
del_make=[5,7, 8, 9, 14,12, 15, 16, 20, 25, 27, 31, 33, 34, 36, 46, 50, 56, 60];
d_make(:,del_make)=[];


%Create market shares
[~, ~, kk] = unique([id_market], 'rows');
out=accumarray(kk,q);
g_mrkt=out(kk);
clear kk out
s0=1-g_mrkt./tm;
sj=q./tm;
deltaold=log(sj)-log(s0); % independent variable

%Create subgroup shares
[~, ~, kk] = unique([id_market, id_fueltype,id_staseg], 'rows');
out=accumarray(kk,q);
g_sgr_t=out(kk);
clear kk out

%Create group shares
[~, ~, kk] = unique([id_market, id_fueltype], 'rows');
out=accumarray(kk,q);
g_gr_f=out(kk);
clear kk out

sjg=q./g_sgr_t;
sgh=g_sgr_t./g_gr_f;
l_sjg=log(sjg);
l_sgh=log(sgh);

% Create explanatory variables and instruments


x_d=[ones(size(p)),kw,consumption,transmission,elec_char_est,drv_rng,turbo, d_make, d_year_month];
%x_d=[ones(size(p)), d_make, d_year_month];

z_d=[table2array(Z),new_inst_seg_sd, x_d];
%z_d(:,[1,3,5,8,11,14,19,20])=[];
z_d(:,[1,3,7,8,9,11,14,15,17,19,20,24])=[];

zz=z_d(:,1:end-size(x_d,2));

Z(:,[1,3,7,8,9,11,14,15,17,19,20,24])=[];

x_lab={'sigma_subg', 'sigma_gr', 'price', 'inter', 'kw', 'consumption', 'transmission', 'electric','char est','Bat','turbo'};
%x_lab={'sigma_subg', 'sigma_gr', 'price'};

% first stage
x_i=[l_sjg l_sgh p x_d];
x_hat = z_d*inv(z_d'*z_d)*z_d'*[l_sjg l_sgh p x_d];
B2sls=inv(x_hat'*x_hat)*x_hat'*deltaold; % nested logit coefficients
mu_hat=[l_sjg l_sgh p x_d]-x_hat;
Var_res=(mu_hat'*mu_hat)/(N-size(B2sls,2));
gamma_par=inv(x_i'*x_i)*x_i'*z_d;
gamma_z=gamma_par(:,1:end-size(x_d,2));
gamma_w=gamma_par(:,end-size(x_d,2)+1:end);
covar_fs=gamma_z*cov(zz)*gamma_z'+gamma_w*cov(x_d)*gamma_w';
covar_x_hat=cov(x_hat);
stderr_gamm  = sqrt(diag(Var_res));

PI=z_d*inv(z_d'*z_d)*z_d';
res_2sls=deltaold-[l_sjg l_sgh p x_d]*B2sls;
OM=(res_2sls'*res_2sls)/length(res_2sls);
covar_2sls=OM*inv([l_sjg l_sgh p x_d]'*PI*[l_sjg l_sgh p x_d]);
sd_2sls=sqrt(diag(covar_2sls)); %Coefficient's stadard deviations.

JJ=size(x_d,1)*(res_2sls'*z_d*inv(z_d'*z_d)*z_d'*res_2sls)/(res_2sls'*res_2sls);

prob_JJ=chi2cdf(JJ,size(z_d,2)-size(x_d,2));

display('The results from the IV estimation: (coeff se)');
resul_IV=table(B2sls(1:size(x_lab,2),:), sd_2sls(1:size(x_lab,2),:), 'RowNames',x_lab)

% Logit
% first stage
x_hatL = z_d*inv(z_d'*z_d)*z_d'*[p x_d];
B2slsL=inv(x_hatL'*x_hatL)*x_hatL'*deltaold; % nested logit coefficients

PIL=z_d*inv(z_d'*z_d)*z_d';
res_2slsL=deltaold-[p x_d]*B2slsL;
OML=(res_2slsL'*res_2slsL)/length(res_2slsL);
covar_2slsL=OML*inv([p x_d]'*PIL*[p x_d]);
sd_2slsL=sqrt(diag(covar_2slsL)); %Coefficient's stadard deviations.

display('The results from the IV estimation: (coeff se)');
resul_IV_L=table(B2slsL(1:size(x_lab,2)-2,:), sd_2slsL(1:size(x_lab,2)-2,:), 'RowNames',x_lab(3:end))
% BLP
% Create supply explanatory matrix

graph_vars=[p,consumption,kw,ccm,bat_size, weight, drv_rng];
descriptive_stats=table((mean(graph_vars)'),(std(graph_vars)'),(prctile(graph_vars,90)'),(prctile(graph_vars,10)'),'VariableNames',{'Mean','Std_dev','Ninetieth_prc','Tenth_prc'}, 'RowNames',{'Price','Fuel_consumption', 'Engine_power','Engine_size', 'Battery_size','Weight','Driving_range'});
writetable(descriptive_stats,'Table_1_Descriptive_stats.txt','Delimiter',' ','WriteRowNames',true) 


x_d=[ones(size(p)),kw,consumption,transmission,elec_char_est,drv_rng,turbo, d_make, d_year_month, d_staseg];
x_s=[ones(size(p)),power,consumption,transmission,bat_size, d_make, d_year_month d_staseg];

z_d=[table2array(Z),new_inst_seg_sd, x_d];
%z_d(:,[1,3,5,8,11,14,19,20])=[];
z_d(:,[1,3,7,8,9,11,14,15,17,19,20])=[];

z_s=[table2array(Z),new_inst_seg_sd,bat_x_lit,tra_x_wei, x_s];
%z_s(:,[1,6,7,11,17,19])=[];
z_s(:,[4,17,19])=[];
% random coeff on: price and fuel consumption
x2=[consumption]; 

k2=size(x2,2)+1; %number of random coefficients
k_S=size(x_s,2);  % number of suply stimation parameters
k_D=size(x_d,2);  % number of suply stimation parameters

% Get random draws
ns=1000;
[v1,vp]=getranddraw(k2-1,(income/1000), ns);
v=[log(vp);v1];

x_lab_d={'price','inter', 'Kw','consumption', 'transmission','Electric','Ele X Char','driving range','turbo'};
        
x_lab_s={'inter', 'power', 'consumption', 'transmission','battery_size'};

% Instruments for supply same as for demand.

%z_s=z_d;

d_d=size(z_d,2);  % number of included and excluded instruments
%z_S=z_D; % supply instruments
d_s=size(z_s,2);  
%getting initial values

ini=[0.5000    0.0001    0.0011];


ztilde=blkdiag(z_d,z_s); 
xtilde=blkdiag(x_d,x_s);
d=size(ztilde,2); 

W=inv(ztilde'*ztilde); % preliminary weight matrix 
%% GMM Demand side estimation 
 
options = optimset('Display','iter-detailed' ,'Algorithm','Interior-Point','GradObj', ...
'off','GradCon','off','DerivativeCheck','off', 'maxiter', 500,...
'MaxFunEvals', 5000);

lb=[0, zeros(1,size(ini,2)-2),0];  %lower bound for mean and random coefficient for price = -inf
ub=[inf, ones(1,size(ini,2)-2)*inf,inf]; %upper bound for mean and random coefficient for price = inf

%[theta2,fval,exitflag,output]=fminsearch(@(thet2) objective_fun(thet2,sj,deltaold,v,x2,p,QI,ind_market,A1,xtilde,ztilde,W,sales_tax),ini,options);

[theta2,fval,exitflag,output]  = fmincon(@(thet)objective_fun(thet,sj,deltaold,v,x2,...
p,QI,ind_market,A1,xtilde,ztilde,W,sales_tax),ini,[],  [], [], [], ...
     [lb],[ub],[], options);

[ya, ma]=computey(theta2,sj,v,x2,p,QI,ind_market,A1,deltaold,sales_tax);
y=[ya;ma];
%W2=inv(z_D'*z_D);
%theta1=inv(x_D'*z_D*W2*z_D'*x_D)*x_D'*z_D*W2*z_D'*ya; % compute linear params theta1=[beta' gamma']'
theta1=inv(xtilde'*ztilde*W*ztilde'*xtilde)*xtilde'*ztilde*W*ztilde'*y; % compute linear params theta1=[beta' gamma']'
resid=y-xtilde*theta1; 
xi=resid(1:N,1);
omega=resid(N+1:end,1);

% Standar errors
Varcovar=covar(theta2,resid,sj,v,x2,p,QI,ind_market,A1,W,d_d,ztilde,xtilde,deltaold,sales_tax);  % compute Var-Covar matrix
se=sqrt(diag(Varcovar));

display('The results for beta are: (coeff se)');
resul_beta=table([theta2(1,1);theta1(1:size(x_lab_d,2)-1,:)],[se(end-k2);se(1:size(x_lab_d,2)-1,:)], 'RowNames',x_lab_d,'VariableNames',{'Ceff','Std_Dev'})
%resul_beta=table(theta1, 'RowNames',x_Dstring(2:end),'VariableNames',{'Ceff'})

display('The results for gamma are: (coeff se)');
resul_gamma=table(theta1(k_D+1:k_D+size(x_lab_s,2),1),[se(k_D+1:k_D+size(x_lab_s,2),1)], 'RowNames',x_lab_s,'VariableNames',{'Ceff','Std_Dev'})

delta=y(1:N,1);
mc_update=exp(ma);
%mc_update=ma;
mark_update = (p./(1+sales_tax))-mc_update;
pcm_update=mark_update./p;
display('The Median markup is:')
median(mark_update)
display('The Median pcm is:')
median(pcm_update)

%%
save('results_2_int')
Elast = elast_blp(theta2,v,x2,p,delta,ind_market,sj,QI);

mean_elast=[];
for jj=1:max(ind_market)
    ttpp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),diag(Elast{:,:,jj}),[],@mean);
    mean_elast=[mean_elast;[repmat(jj,size(ttpp)),ttpp]];     
end
M_E=reshape(mean_elast(:,2),4,size(mean_elast,1)/4); 
Mean_Elast = sum(M_E,2) ./ sum(M_E~=0,2)


% Marginal effect of power and power

Marg_con = marginal_effect(theta2,theta1(3),theta2(2),v(2,:),v,x2,consumption,delta,ind_market,sj,QI);

mean_marg_con=[];
for jj=1:max(ind_market)
    ttpp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),diag(Marg_con{:,:,jj}),[],@mean);
    mean_marg_con=[mean_marg_con;[repmat(jj,size(ttpp)),ttpp]];     
end
M_M_C=reshape(mean_marg_con(:,2),4,size(mean_marg_con,1)/4); 
Mean_ME_C = sum(M_M_C,2) ./ sum(M_M_C~=0,2)


Marg_pow = marginal_effect(theta2,theta1(2),0,v(2,:),v,x2,kw,delta,ind_market,sj,QI);

mean_marg_pow=[];
for jj=1:max(ind_market)
    ttpp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),diag(Marg_pow{:,:,jj}),[],@mean);
    mean_marg_pow=[mean_marg_pow;[repmat(jj,size(ttpp)),ttpp]];     
end
M_M_P=reshape(mean_marg_pow(:,2),4,size(mean_marg_pow,1)/4); 
Mean_ME_P = sum(M_M_P,2) ./ sum(M_M_P~=0,2)

% Willingness to pay consumption
Marg_price = marginal_effect(theta2,-theta2(1),theta2(end),v(1,:),v,x2,p,delta,ind_market,sj,QI);

mean_wtp=[];
for jj=1:max(ind_market)
    ttpp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),diag(Marg_con{:,:,jj})./(diag(Marg_price{:,:,jj})),[],@mean);
    mean_wtp=[mean_wtp;[repmat(jj,size(ttpp)),ttpp]];     
end
M_WP=reshape( mean_wtp(:,2),4,size( mean_wtp,1)/4); 
Mean_WTP= sum(M_WP,2) ./ sum(M_WP~=0,2)

wtp_power=(Mean_ME_P./Mean_WTP*1000)
wtp_gasoline=(Mean_ME_C/Mean_WTP*1000)


%% Check elasticities
CPE_E=[];
for jj=1:max(ind_market)
    [old_sj_fuel,new_sj_fuel,sj_change_fuel] = change_ms(sj,p,Elast,ind_market,id_fuelcat,electric,0.1, jj);
    CPE_E=[CPE_E;sj_change_fuel];
end
CPE_E(isnan(CPE_E))=0;
CP_EE=reshape(CPE_E,4,size(CPE_E,1)/4);
Mean_CPE_E = sum(CP_EE,2) ./ sum(CP_EE~=0,2);

[~,~,mf]=unique([id_market,id_fuelcat], 'rows');

ms_ini=accumarray(mf,sj);
p_new=p.*(1+0.1.*(electric));

expmu_new=expmu(theta2,v,x2,p_new);      
sj_new=mktshares(delta,expmu_new,QI,ind_market);

ms_fin=accumarray(mf,sj_new);

change=[ms_ini, ms_fin];


%% Estimation with updatd weight matrix

ini=theta2; 
clear theta2;
TT=[repmat(xi,1,d_d).*z_d repmat(omega,1,d_s).*z_s];
V=TT'*TT;
W_update=inv(V); 

[theta2_updated,fval_updated,exitflag,output]  = fmincon(@(thet)objective_fun(thet,sj,deltaold,v,x2,...
p,QI,ind_market,A1,xtilde,ztilde,W_update,sales_tax),ini,[],  [], [], [], ...
     [lb],[ub],[], options);

[ya, ma]=computey(theta2_updated,sj,v,x2,p,QI,ind_market,A1,deltaold,sales_tax);
y=[ya;ma];
theta1_updated=inv(xtilde'*ztilde*W_update*ztilde'*xtilde)*xtilde'*ztilde*W_update*ztilde'*y; % compute linear params theta1=[beta' gamma']'
resid_updated=y-xtilde*theta1_updated; 
xi_updated=resid_updated(1:N,1);
omega_updated=resid_updated(N+1:end,1);

Varcovar_update=covar(theta2_updated,resid_updated,sj,v,x2,p,QI,ind_market,A1,W,d_d,ztilde,xtilde,deltaold,sales_tax);  % compute Var-Covar matrix
se_update=sqrt(diag(Varcovar_update));

display('The results for beta are: (coeff se)');
resul_beta=table([theta2_updated(1,1);theta1_updated(1:size(x_lab_d,2)-1,:)],[se_update(end-k2);se_update(1:size(x_lab_d,2)-1,:)], 'RowNames',x_lab_d,'VariableNames',{'Ceff','Std_Dev'})
%resul_beta=table(theta1, 'RowNames',x_Dstring(2:end),'VariableNames',{'Ceff'})

display('The results for gamma are: (coeff se)');
resul_gamma=table(theta1_updated(k_D+1:k_D+size(x_lab_s,2),1),[se_update(k_D+1:k_D+size(x_lab_s,2),1)], 'RowNames',x_lab_s,'VariableNames',{'Ceff','Std_Dev'})

writetable(resul_beta,'Random_Coef_demand.txt','Delimiter',' ','WriteRowNames',true)  
writetable(resul_gamma,'Random_Coef_supply.txt','Delimiter',' ','WriteRowNames',true)  

delta=y(1:N,1);
mc_update=exp(ma);
%mc_update=ma;
mark_update = (p./(1+sales_tax))-mc_update;
pcm_update=mark_update./p;
display('The Median markup is:')
median(mark_update)
display('The Median pcm is:')
median(pcm_update)
save('results_new_instruments')



expmu1=expmu(theta2_updated,v,x2,p);
sj_hat=mktshares(delta,expmu1,QI,ind_market);
RSS_sj=(sj-sj_hat)'*(sj-sj_hat);
TSS_sj=(sj-mean(sj))'*(sj-mean(sj_hat));

RSS_delta=(deltaold-delta)'*(deltaold-delta);
TSS_delta=(deltaold-mean(deltaold))'*(deltaold-mean(deltaold));

R_sq_sj=1-(RSS_sj/TSS_sj)
MSE_sj=RSS_sj./length(sj)
MSE_delta=RSS_delta./length(delta)

%BPT

aux_xi = fitlm(x_d(:,2:end),xi_updated.^2);
df_xi = size(x_d(:,2:end),2); 
n_test = numel(sj(2:end));
T_xi = aux_xi.Rsquared.Ordinary*n_test;
P_xi = 1-chi2cdf(abs(T_xi),df_xi);

aux_omega = fitlm(x_s(:,2:end),omega_updated.^2);
df_omega = size(x_s(:,2:end),2); 
T_omega = aux_omega.Rsquared.Ordinary*n_test;
P_omega = 1-chi2cdf(abs(T_omega),df_omega)


histogram(xi_updated)
histogram(omega_updated)
%% Get elasticity

Elast = elast_blp(theta2_updated,v,x2,p,delta,ind_market,sj,QI);

mean_elast=[];
for jj=1:max(ind_market)
    ttpp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),diag(Elast{:,:,jj}),[],@mean);
    mean_elast=[mean_elast;[repmat(jj,size(ttpp)),ttpp]];     
end
M_E=reshape(mean_elast(:,2),4,size(mean_elast,1)/4); 
Mean_Elast = sum(M_E,2) ./ sum(M_E~=0,2)

Mean_prof_mar=accumarray(grp2idx(id_fuelcat),pcm_update,[],@mean);
Mean_markup=accumarray(grp2idx(id_fuelcat),mark_update,[],@mean);

save('results_final')

% Elast per consumption cuant
mean_elast_cons=[];
for jj=1:max(ind_market)
    ttpp_c=accumarray(grp2idx(cons_cat(ind_market==jj)),diag(Elast{:,:,jj}),[],@mean);
    ttpp_d=accumarray(grp2idx(cons_cat(ind_market==jj)),cons_cat(ind_market==jj),[],@mean);
    mean_elast_cons=[mean_elast_cons;[repmat(jj,size(ttpp_c)),ttpp_c,ttpp_d]];     
end

[~, ~, ii ] = unique(mean_elast_cons(:,3));
Mean_Elast_cons = accumarray(ii, mean_elast_cons(:,2), [], @mean);
Median_Elast_cons = accumarray(ii, mean_elast_cons(:,2), [], @median);

% Markup by price quantile
mean_mk_pquant=[];
for jj=1:max(ind_market)
    ttppp_c=accumarray(grp2idx(pr_cat(ind_market==jj)),mark_update(ind_market==jj),[],@mean);
    ttppp_d=accumarray(grp2idx(pr_cat(ind_market==jj)),pr_cat(ind_market==jj),[],@mean);
    mean_mk_pquant=[mean_mk_pquant;[repmat(jj,size(ttppp_c)),ttppp_c,ttppp_d]];     
end

[~, ~, ii ] = unique(mean_mk_pquant(:,3));
Mean_mkup_p = accumarray(ii, mean_mk_pquant(:,2), [], @mean);
Median_mkup_p = accumarray(ii, mean_mk_pquant(:,2), [], @median);

Elast_makup=table(Mean_Elast,Mean_prof_mar,Mean_markup*1000, 'RowNames',{'Electric','Natural Gas','Hybrid','Gasoline'},'VariableNames',{'Mean_own_price_elasticity','Mean_profit_margin','Mean_markup_USD'});
writetable(Elast_makup,'Elast_Mark_up.txt','Delimiter',' ','WriteRowNames',true)  

% Market share by vehicle type
 ms_elec=accumarray(ind_market,sj.*electric)./ms_markt;
 ms_hybrid=accumarray(ind_market,sj.*hybrid)./ms_markt;
 year_mon=accumarray(ind_market,year,[],@median);
 
 ms_elec_y=accumarray(year_mon,ms_elec,[],@median);
 ms_hybrid_y=accumarray(year_mon,ms_hybrid,[],@median);
 
small_ms=sj.*small_car;
mean_ms_small=accumarray(ind_ver,small_ms,[],@mean);
[ii,ii] = sort(mean_ms_small);
sml_1=unique(id_ver(ind_ver==ii(end)))
mean(sj(id_ver==sml_1))./mean(ms_markt)
sml_2=unique(id_ver(ind_ver==ii(end-1)))
mean(sj(id_ver==sml_2))./mean(ms_markt)

mini_car=id_staseg=='1 Mini';
mini_ms=sj.*mini_car;
mean_ms_mini=accumarray(ind_ver,mini_ms,[],@mean);
[ii,ii] = sort(mean_ms_mini);
mini_1=unique(id_ver(ind_ver==ii(end)))
mean(sj(id_ver==mini_1))./mean(ms_markt)
mini_2=unique(id_ver(ind_ver==ii(end-1)))
mean(sj(id_ver==mini_2))./mean(ms_markt)

%% Price elasticity by gas types

expmu_ori=expmu(theta2_updated,v,x2,p);
sj_ori=mktshares(delta,expmu_ori,QI,ind_market);

p2_gas=p.*(ones(size(p))+(0.10*(id_fuelcat=='NR')));
expmu_gas=expmu(theta2_updated,v,x2,p2_gas);
sj_gas=mktshares(delta,expmu_gas,QI,ind_market);

change_md_gas=(sj_gas./sj)-1;
Mean_chang_gas=accumarray(grp2idx(id_fuelcat),change_md_gas,[],@mean);
Mean_elas_gas=Mean_chang_gas/.1;

p2_hyb=p.*(ones(size(p))+(0.10*(id_fuelcat=='Hybrid')));
expmu_hyb=expmu(theta2_updated,v,x2,p2_hyb);
sj_hyb=mktshares(delta,expmu_hyb,QI,ind_market);

change_md_hyb=(sj_hyb./sj)-1;
Mean_chang_hyb=accumarray(grp2idx(id_fuelcat),change_md_hyb,[],@mean);
Mean_elas_hyb=Mean_chang_hyb/.1;

p2_ele=p.*(ones(size(p))+(0.10*(id_fuelcat=='Electric')));
expmu_ele=expmu(theta2_updated,v,x2,p2_ele);
sj_ele=mktshares(delta,expmu_ele,QI,ind_market);

change_md_ele=(sj_ele./sj)-1;
Mean_chang_ele=accumarray(grp2idx(id_fuelcat),change_md_ele,[],@mean);
Mean_elas_ele=Mean_chang_ele/.1;

Elast_gas_type=table(Mean_elas_ele,Mean_elas_hyb,Mean_elas_gas, 'RowNames',{'Electric','Natural Gas','Hybrid','Gasoline'},'VariableNames',{'Electric','Hybrid','Gasoline'});
writetable(Elast_gas_type,'Elast_Gas_type.txt','Delimiter',' ','WriteRowNames',true)  


%% Get owm and cross price elasticities

sigma_sg=B2sls(1);
sigma_g=B2sls(2);
alpha=B2sls(3);

delta=[p x_d]*B2sls(3:end);

parameters=[sigma_g,sigma_sg];

grouping=[id_market,id_fueltype,id_staseg];

S=marketshare(parameters,delta,grouping);

[elst, mrkp]=elast(S(:,1),S(:,2),S(:,3),alpha,sigma_g,sigma_sg,p,grouping,A1); % Compute own and cross price elasticities.

mc=(p-mrkp)./(1+tariffs); % recover marginal cost

% Compute cross elasticities

N=size(p);
temp=repmat(grp2idx(id_market),[1,N]);
MK=eq(temp,temp');


% Inside groups

mean_elst_el_el=sum(sum(electric.*elst.*electric'.*(1-eye(length(p)))))./sum(sum((electric.*ones(size(elst)).*electric'.*(1-eye(length(p)))).*MK))
mean_elst_hy_hy=sum(sum(hybrid.*elst.*hybrid'.*(1-eye(length(p)))))./sum(sum((hybrid.*ones(size(elst)).*hybrid'.*(1-eye(length(p)))).*MK))
mean_elst_nr_nr=sum(sum(non_elec.*elst.*non_elec'.*(1-eye(length(p)))))./sum(sum((non_elec.*ones(size(elst)).*non_elec'.*(1-eye(length(p)))).*MK))


% Elect

mean_elst_el_hy=sum(sum(electric.*elst.*hybrid'.*(1-eye(length(p)))))./sum(sum((electric.*ones(size(elst)).*hybrid'.*(1-eye(length(p)))).*MK))
mean_elst_el_nr=sum(sum(electric.*elst.*non_elec'.*(1-eye(length(p)))))./sum(sum((electric.*ones(size(elst)).*non_elec'.*(1-eye(length(p)))).*MK))

% Hybrid

mean_elst_hy_el=sum(sum(hybrid.*elst.*electric'.*(1-eye(length(p)))))./sum(sum((hybrid.*ones(size(elst)).*electric'.*(1-eye(length(p)))).*MK))
mean_elst_hy_nr=sum(sum(hybrid.*elst.*non_elec'.*(1-eye(length(p)))))./sum(sum((hybrid.*ones(size(elst)).*non_elec'.*(1-eye(length(p)))).*MK))

% NR

mean_elst_nr_el=sum(sum(non_elec.*elst.*electric'.*(1-eye(length(p)))))./sum(sum((non_elec.*ones(size(elst)).*electric'.*(1-eye(length(p)))).*MK))
mean_elst_nr_hy=sum(sum(non_elec.*elst.*hybrid'.*(1-eye(length(p)))))./sum(sum((non_elec.*ones(size(elst)).*hybrid'.*(1-eye(length(p)))).*MK))
% Base line Market share

[gr2, ~, jj] = unique([id_market,id_fuelcat] , 'rows');
year_2=accumarray(jj,year,[],@mean);

ms_bl=accumarray(jj,prod(S,2));
[gr3, ~, kk] = unique([gr2(:,2),categorical(year_2)] , 'rows');

ms_bl2=accumarray(kk,ms_bl,[],@mean);

% Mean own price elast

own_pe=diag(elst);

mean_own_pe=accumarray(grp2idx(id_fuelcat),own_pe,[],@mean);

%% Simulation
tariffs_cf=ones(size(p))*.35;

p_cf=mc.*(1+tariffs_cf)+mrkp; 

delta_cf=[p_cf x_d]*B2sls(3:end);

S_cf=marketshare(parameters,delta_cf,grouping);

[elst_cf, mrkp_cf]=elast(S_cf(:,1),S_cf(:,2),S_cf(:,3),alpha,sigma_g,sigma_sg,p_cf,grouping,A1);

ms_cf=accumarray(jj,prod(S_cf,2));
ms_cf2=accumarray(kk,ms_cf,[],@mean);

chang_ms=1-ms_cf2./ms_bl2;

% Consumer's and Producer's surplus

[~, ~, gg] = unique(id_market , 'rows');
exdelta_bl=accumarray(gg,exp(delta));
exdelta_cf=accumarray(gg,exp(delta_cf));


CS_bl=mean(((0.577+(log(1+exdelta_bl)))./-alpha),2);
CS_cf=mean(((0.577+(log(1+exdelta_cf)))./-alpha),2);

chang_CS=(1-CS_cf./CS_bl);

% Producer's and Producer's surplus

margin_bl=mrkp;
margin_cf=mrkp_cf;

PS_bl=accumarray(gg,margin_bl);
PS_cf=accumarray(gg,margin_cf);

chang_PS=1-(PS_cf./PS_bl);

% Calculation emissions

ave_millage=20000./1.60934;

Units_bl=prod(S,2).*(1-s0).*tm;
Units_cf=prod(S_cf,2).*(1-s0).*tm;

emi_bl=emissions.*Units_bl.*ave_millage./100;
emi_cf=emissions.*Units_cf.*ave_millage./100;

T_emi_bl=accumarray(gg,emi_bl);
T_emi_cf=accumarray(gg,emi_cf);

chang_EM=1-(T_emi_cf./T_emi_bl);

% Calculate tax revenue

tax_bl=mc.*tariffs.*Units_bl;
tax_cf=mc.*tariffs_cf.*Units_cf;

T_tax_bl=accumarray(gg,tax_bl);
T_tax_cf=accumarray(gg,tax_cf);

chang_tax=1-(T_tax_cf./T_tax_bl);
